In [1]:
import coldatoms
import numpy as np
%matplotlib notebook
import matplotlib.pyplot as plt

Introduction

Doppler cooling is one of the most important experimental techniques in cold atom science. Perhaps it's indicative of the impact of this technology that at least five of the names associated with the early days of laser cooling (Wineland, Hansch, Chu, Tannoudji, and Phillips) went on to earn Nobel Prizes in Physics, most of them directly for laser cooling and others for the lifetime contributions to cold atom science. William D. Phillips' Nobel lecture gives a good overview of some of the ways in which laser cooling has had an impact in low energy physics and other areas of physics research.

In this notebook we explore how Doppler cooling, the simplest form of laser cooling, can be modeled in the coldatoms library.

Central to laser cooling is the radiation pressure force generated by photons resonantly scattering off of atoms. For a two level atom the scattering rate is given by

$$ \dot{n} = S\frac{\gamma_0}{2\pi}\frac{\left(\gamma_0/2\right)^2}{(\gamma_0/2)^2(1+2S)+\Delta^2(\mathbf{v})} $$

In this equation, $S$ is the intensity of the laser in units of the saturation intesnity, $2\pi/\gamma_0$ is the excited state lifetime, and $\Delta$ is the detuning of the laser frequency from the atomic resonance frequency.

Each time the atom absorbs and reemits a photon it receives a recoil kick. Very absorbed photon comes out the the beam and therefore has a well defined momentum ($\hbar \mathbf{k}$ where $\mathbf{k}$ is the wavevector of the laser). The emitted photons travel in more or less random directions and the force due to them therefore averages approximately to zero. The net result is a forc eon the atom along the direction of propagagation of the laser beam and a fluctuating force that is more or less isotropic.

Now comes the imporatnt part that allows us to use this force for cooling: The detuning $\Delta$ of the laser from the atomic transition is velocity dependent due to the Doppler shift. In free space we have

$$ \Delta(\mathbf{v}) = \omega_L - \omega_A - \mathbf{k}\cdot\mathbf{v} $$

If we then "red detune" the laser i.e. we choose a laser frequency such that $\omega_L<\omega_A$, it is easy to see that atoms moving towards the laser with $\mathbf{k}\cdot\mathbf{v}<0$, will scatter more photons than atoms moving away from the laser. They will hence experience a decelerating force if they're moving towards the laser beam.

"Wait a second" you say. "That's all good and fine. The atoms are slowed down if they're moving in the opposite direction as the laser's propagation direction. But what if they move in the direction of propagation of the laser. Won't they get accelerated in that case?"

You are correct! One laser beam can only slow down the atoms if they're going one way. To slow them down going either direction we need a second laser that is propagating in the opposite direction. By combining six such lasers, a pair for each Cartesian direction, we can achieve cooling of the atoms' motion in all three spatial directions.

We have neglected the fluctuating force due to the emitted photons so far. Unfortunately these recoil kicks never cancel completely because they are uncorrelated with one another. The recoil kicks make th atoms undergo a random walk in momentum space. They are a heating mechanism.

The balance between the cooling rate due to the coeherent friction force from absorption and the heating due to the incoherent emission events corresponds to the lowest temperature that can be achieved by Doppler cooling. This temperature is called the Doppler temperature. We will determine it computationally later in this notebook.

It may be worth mentioning that many laser cooling schemes exist that are able to achieve temperatures lower than the Doppler limit. We will not consider these so called sub-Doppler schemes here.

Doppler cooling in coldatoms

So how dow we simulate Doppler cooling using the coldatoms library? The answer is we use the RadiationPressure force. This force mimicks the momentum recoil due to absorption and emission of photons in resonance fluorescence. The RadiationPressure force is completely determined by the excited state decay rate $\gamma_0$, the driving laser's wavevector $\mathbf{k}$, the laser intensity, and the detuning of the laser from the atomic transition frequency. The intensity is a function of the atomic position because the laser intensity varies throughout space. The detuning depends on the atomic position and velocity. It depends on position because external fields may lead to shifts of the atomic transition frequency (e.g. magnetic fields leading to Zeeman shifts) and it depends on velocity via the Doppler shift.

In this notebook we consider a well collimated Gaussian laser beam. It has an intensity profile given by

$$ S(\mathbf{x})=S_0e^{-x_\perp^2/\sigma^2} $$

where $S_0$ is the peak intensity of the beam, $\sigma$ is the width of the beam, and $x_\perp=\mathbf{x}-\mathbf{x_0} - (\mathbf{x}-\mathbf{x_0})\cdot \mathbf{k}/k$ is the distance of the atom from the center of the beam. We represent the intensity by the following Python class:


In [2]:
class GaussianBeam(object):
    """A laser beam with a Gaussian intensity profile."""
    
    def __init__(self, S0, x0, k, sigma):
        """Construct a Gaussian laser beam from position, direction, and width.
        
        S0 -- Peak intensity (in units of the saturation intensity).
        x0 -- A location on the center of the beam.
        k -- Propagation direction of the beam (need not be normalized).
        sigma -- 1/e width of the beam."""
        self.S0 = S0
        self.x0 = np.copy(x0)
        self.k_hat = k / np.linalg.norm(k)
        self.sigma = sigma
        
    def intensities(self, x):
        xp = x - self.x0
        xperp = xp - np.outer(xp.dot(self.k_hat[:, np.newaxis]), self.k_hat)
        return self.S0 * np.exp(-np.linalg.norm(xperp, axis=1)**2/self.sigma**2)

The following figure shows a contour plot of such a beam originating at $\mathbf{x}_0=(1,1,0)^T$ and propagating in the $\mathbf{k}=(1, 1, 0)^T$ direction.


In [3]:
beam = GaussianBeam(1.0, np.array([1.0,1.0,0.0]),np.array([1.0, 1.0, 0.0]), 1.0)

num_pts = 10
x_min = -3
x_max = 3
y_min = -3
y_max = 3

x = np.linspace(x_min, x_max, num_pts)
y = np.linspace(x_min, x_max, num_pts)
pts = np.array([[x[i], y[j], 0] for i in range(num_pts) for j in range(num_pts)])
intensities = beam.intensities(pts).reshape(num_pts, num_pts)
xx, yy = np.meshgrid(x, y)
plt.contour(xx, yy, intensities)
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')


Out[3]:
<matplotlib.text.Text at 0x7f74391b5da0>

Besides the intensity we also need to tell RadiationPressure what the laser-atom detuning is. Here we assume that we only need to account for Doppler shifts and thus we have

$$ \Delta(\mathbf{x}, \mathbf{v}) = \Delta_0-\mathbf{k}\cdot\mathbf{v}. $$

The frequency $\Delta_0$ is the detuning between laser and atomic transition when the atom is at rest. Here is a class to represent this detuning:


In [4]:
class DopplerDetuning(object):
    def __init__(self, Delta0, k):
        self.Delta0 = Delta0
        self.k = np.copy(k)
        
    def detunings(self, x, v):
        return self.Delta0 - np.inner(self.k, v)

In [5]:
detuning = DopplerDetuning(-0.5, (1, 0, 0))

One dimensional laser cooling

As a first example we consider a single atom being laser cooled along the $x$ dimension. We take ${}^{87}\rm{Rb}$:


In [53]:
ensemble = coldatoms.Ensemble(num_ptcls=1)
ensemble.ensemble_properties['mass'] = 87*1.67e-27
wavelength = 780.0e-9
k = 2.0 * np.pi / wavelength
gamma = 2.0 * np.pi * 6.1e6
hbar = 1.0e-34
sigma = 1.0e-3

To represent the 1D MOT we need to add two radiation pressure forces. One for the laser propagating along the $+x$ direction and one propagating along the $-x$ direction. We use a an intensity of $S=0.1S_0$ and a beam width of $\sigma = 1{\rm mm}$:


In [54]:
one_d_mot = [
    coldatoms.RadiationPressure(gamma, np.array([hbar * k, 0.0, 0.0]),
                      GaussianBeam(S0=0.1, x0=np.array([0.0, 0.0, 0.0]), k=np.array([k, 0.0, 0.0]), sigma=sigma),
                      DopplerDetuning(-0.5 * gamma, np.array([k, 0.0, 0.0]))),
    coldatoms.RadiationPressure(gamma, np.array([-hbar * k, 0.0, 0.0]),
                      GaussianBeam(S0=0.1, x0=np.array([0.0, 0.0, 0.0]), k=np.array([-k, 0.0, 0.0]), sigma=sigma),
                      DopplerDetuning(-0.5 * gamma, np.array([-k, 0.0, 0.0])))]

We have the atom start with a velocity of $v_x = 5\rm{m}/\rm{s}$ and we simulate its evolution with three different time step sizes to check if our simulation is converged:


In [55]:
velocities = []
times = []

# Loop over time step sizes.
for i in range(3):
    # Reset positions and velocities to initial conditions.
    ensemble.x *= 0.0
    ensemble.v *= 0.0
    ensemble.v[0, 0] = 5.0e0
    
    # The time step size.
    dt = 1.0e-5 / 2**i
    
    v = []
    t = []
    # Now do the time integration and record velocities and times.
    for i in range(2000 * 2**i):
        v.append(ensemble.v[0, 0])
        t.append(i * dt)
        coldatoms.drift_kick(dt=dt, ensemble=ensemble, forces=one_d_mot)
    velocities.append(v)
    times.append(t)

The following plot shows the evolution of the particle's velocity.


In [42]:
plt.figure()
for i in range(3):
    plt.plot(1.0e3 * np.array(times[i]), velocities[i])
plt.xlabel(r'$t/\rm{ms}$')
plt.ylabel(r'$v_x/(\rm{m}/\rm{s})$')


Out[42]:
<matplotlib.text.Text at 0x7f74384c9470>

After about $5\rm{ms}$ the particle has been brought to an almost complete stop. The non-exponential shape of the velocity evolution is due to the finite capture range of the cooling lasers. When the atomic velocity is too large the laser is too far detuned from the atomic transition. The atom then doesn't scatter any photons and hence does not get decelerated.

The three simulation traces agree well with one another indicating that the numerical integration is converged.

3D optical molasses

Next we consider fully three dimensional laser cooling as is often used in magneto-optical traps. To obtain cooling in all three dimensions we need three pairs of lasers:


In [44]:
three_d_mot = [
    coldatoms.RadiationPressure(gamma, hbar * kp,
                      GaussianBeam(S0=0.1, x0=np.array([0.0, 0.0, 0.0]), k=kp, sigma=1.0e-3),
                      DopplerDetuning(-0.5 * gamma, kp)) for kp in [
        np.array([k, 0.0, 0.0]),
        np.array([-k, 0.0, 0.0]),
        np.array([0.0, k, 0.0]),
        np.array([0.0, -k, 0.0]),
        np.array([0.0, 0.0, k]),
        np.array([0.0, 0.0, -k]),
    ]]

The integration of the atoms' equations of motion is virtually unchanged. We merely have to use the forces due to the three dimensional mot:


In [82]:
velocities_3d = []
times_3d = []

# Loop over time step sizes.
for i in range(3):
    # Reset positions and velocities to initial conditions.
    ensemble.x *= 0.0
    ensemble.v *= 0.0
    ensemble.v[0, 0] = 5.0e0
    
    # The time step size.
    dt = 1.0e-5 / 2**i
    
    v = []
    t = []
    # Now do the time integration and record velocities and times.
    for i in range(3000 * 2**i):
        v.append(ensemble.v[0, 0])
        t.append(i * dt)
        coldatoms.drift_kick(dt=dt, ensemble=ensemble, forces=three_d_mot)
    velocities_3d.append(v)
    times_3d.append(t)

When we now consider the $x$ component of the atomic velocity we see that the residual velocity fluctuations are larger. This is because now noise from two additional pairs of lasers feed into the $x$ component:


In [83]:
plt.figure()
for i in range(3):
    plt.plot(1.0e3 * np.array(times_3d[i]), velocities_3d[i])
plt.xlabel(r'$t/\rm{ms}$')
plt.ylabel(r'$v_x/(\rm{m}/\rm{s})$')


Out[83]:
<matplotlib.text.Text at 0x7f74380265c0>

The residual velocity fluctuations correspond to the so-called Doppler temperature. They are caused by the atoms absorbing a photon from the "wrong" beam. We find for the standard deviation of the $x$-component of the velocity:


In [85]:
tmin = [500, 1000, 2000]
for i in range(3):
    print(np.std(np.array(velocities_3d[i][tmin[i]:])))


0.0735762179891
0.0728011230558
0.0794315210151

We can compare that with the value that theory predicts based on the Doppler temperature $\sqrt{\langle v_x^2\rangle}$:


In [60]:
np.sqrt(hbar * gamma / ensemble.ensemble_properties['mass'] / 2.0 / 3.0)


Out[60]:
0.06630730314379045

I think our atoms are slightly too hot but the discrepancy between theory and computation could be due to sampling errors. I'll have to think about other possible explanations. Hope it's not an issue with the random number generator...


In [ ]: